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' The optimal control of problems that are constrained by partial differential equations with un- 

C^ certainties and with uncertain controls is addressed. The Lagrangian that defines the problem 

vN is postulated in terms of stochastic functions, with the control function possibly decomposed 

^_^ into an unknown deterministic component and a known zero-mean stochastic component. The 

extra freedom provided by the stochastic dimension in defining cost functionals is explored, 

demonstrating the scope for controlling statistical aspects of the system response. One-shot 

stochastic finite element methods are used to find approximate solutions to control problems. 

(-H It is shown that applying the stochastic collocation finite element method to the formulated 

- *— * problem leads to a coupling between stochastic collocation points when a deterministic optimal 

JI^ control is considered or when moments are included in the cost functional, thereby forgoing 

\^ the primary advantage of the collocation method over the stochastic Galerkin method for 

the considered problem. The application of the presented methods is demonstrated through 

-^^ a number of numerical examples. The presented framework is sufficiently general to also 

►^ consider a class of inverse problems, and numerical examples of this type are also presented. 

N Keywords: Optimal control, uncertainty, stochastic finite element method, stochastic inverse 

Qv problems, stochastic partial differential equations. 
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t^ 1 Introduction 

o 

'"^ In many applications, forces or boundary conditions are to be determined such that the response of 

, , a physical or engineering system is optimal in some sense. These problems can often be formulated 

^ as the minimisation of an objective functional subject to a set of constraint equations in the form 

k> of partial differential equations (PDEs). For problems that involve uncertainty, incorporating 

V^ stochastic information into a control formulation can lead to a quantification of the statistics of 

C^ the system response. Moreover, there is scope to control a system not only for an optimal mean 

response, but also to include statistics of the response in a cost functional. 'We note that stochastic 

PDE-constrained optimisation problems are closely related to stochastic inverse problems, where 

the control variable corresponds to the parameter to be identified [28l [131 El] • 

We examine in this work the numerical solution of optimal control problems constrained by 
stochastic PDEs and with uncertain controls. Stochastic finite element-based solvers have been 
studied extensively for a large range of stochastic PDEs j24j. However, few results and examples 
on solving optimisation problems constrained by stochastic PDEs are available. Existence of an 
optimal solution to stochastic optimal control problems constrained by stochastic elliptic PDEs 
was studied by Hou et al. [Hj for deterministic control functions. Borzi and von Winckel [7] 
and Borzi studied multigrid solvers for stochastic collocation solutions of parabolic and elliptic 
optimal control problems with random coefficients and stochastic control functions. We contend 



that problems with an unknown stochastic 'control function' constitute stochastic inverse problems 
and are different from control problems where the focus is on computing a deterministic component 
of the control function which forms the control 'signal'. If the stochastic properties of the control 
are computed, ad hoc procedures are required to extract a deterministic function, which will in 
general not be the optimal control. We choose to split the control function into an unknown 
deterministic component, which is to be computed, and a known zero-mean stochastic part that 
represents the uncertainty in the controller response. 

Control cost functionals will be formulated in terms of norms that include both spatial and 
stochastic dimensions. The inclusion of the stochastic dimension provides additional freedom in the 
definition of cost functionals. We formulate a one-shot approach to control, and solve the resulting 
equations via a stochastic collocation or a Galerkin finite element method. The one-shot approach 
is in contrast with methods that require iteration |10| . Overviews of the stochastic collocation 
and Galerkin finite element methods can be found in 3J and [I], respectively. The stochastic 
collocation method is often preferred over the Galerkin approach as it converts a stochastic problem 
into a collection of decoupled deterministic problems. However, we will show that the so-called 
non-intrusivity property of the collocation method that is often exploited, including for stochastic 
inverse problems [28], \7] , does not hold for a large class of stochastic PDE-constrained optimisation 
problems. The stochastic Galerkin method, on the other hand, can be applied straightforwardly to 
stochastic optimal control problems. The efficient solution of stochastic finite element problems, 
and in particular for the stochastic Galerkin method, can hinge on the development and application 
of effective preconditioncrs. This aspect is addressed for stochastic control problems, with two 
preconditioners that take the specific structure of the Galerkin one-shot systems into account 
presented. Extensive numerical examples support our investigations, and the computer code to 
reproduce all numerical results is available under the GNU Lesser Public License (LGPL) as part 
of the supporting material [50] . 

The remainder of this paper is organised as follows. Section [2] presents a control problem 
involving an elliptic stochastic PDE constraint and formulates this problem as a coupled system 
of stochastic PDEs. The framework developed in Section[2]is formulated such that it is sufficiently 
general to also address a class of inverse problems. The stochastic finite element discretisation is 
presented in Section [31 Section [4] considers additional regularising terms and their impact in the 
context of stochastic finite element solvers. Iterative solvers and preconditioners for the one-shot 
Galerkin system are discussed in Section [s] which is followed in Section [6] by numerical examples 
of stochastic optimal control problems. Numerical examples illustrating the solution of stochastic 
inverse problems are given in Section [7| and conclusions are drawn in Section [8] 

2 A control problem with stochastic PDE constraints 

We consider optimal control problems constrained by partial differential equations with stochastic 
coefficients. In general, these can be formulated as: 

min J^{z,u) subject to c(z,u)~0, (1) 

where ^JiZxC/— s-Misa cost (or objective) functional, c is a constraint, z the state variable, u 
the control variable, and Z and U are suitably defined function spaces. 

In this work, the constraint equations in M will be given by one or more PDEs with stochastic 
coefficients. The state variable will be a stochastic function and the control can be either deter- 
ministic or stochastic. Although the cost functional J^ contains stochastic functions, it will be 
defined such that its outcome is deterministic. Before presenting concrete examples, it is useful to 
provide some definitions. 

2.1 Definitions 

We will consider problems posed on a polygonal spatial domain D CW^ with boundary dD and 
where 1 <d <?>. The boundary is decomposed into dDo and SD^r such that ODd U dDif = dD 



and dDo H dD^ = 0. The outward unit normal vector to dD is denoted by n. Consider also 
a complete probability space {Q,J^,V), where fi is a sample space, J^ is a cr-algebra and P is a 
probability measure. We define the tensor product Hilbert space L^ {D) (g) L^ (fi) of second-order 
random fields as 

L^{D)®L^{n)^ [z:D®n^R, f f \z\^dxdV <oo|. (2) 

This space is equipped with the norm 
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uen, 


(6) 



l^llL^(D)«L^(a)=(^y^J^kl'dxd7'j . (3) 

Analogously, the tensor product spaces Hq{D) (E) L'^(ri) and H^{D) (E) i^(51) can be defined [T]. 

2.2 Model problem 

2.2.1 Constraint equation 

As a model constraint, we consider a stochastic steady-state diffusion equation; 

— Va; • (k {x, w) Va;Z {x, a;)) — u (x, w) 
z (x, Lo) = Z]j (a;, uj) 
K (x, w) VxZ {x,uj) ■ n = g (x, lu) 

where the function z : D x fi -h- M is to be found, the diffusion coefhcient n : D y. Q, ^^ M. \s 
prescribed, u : D x Vt ^f M. \s a. source term, zo : dDo x fi — >■ M is a Dirichlet condition and 
g : dD^ x 57 ^ M is a Neumann condition. The operator V^: involves only derivatives with 
respect to the spatial variable x. If k, m, zd and g are second-order random fields and assuming 
boundedness and positivity of the diffusion coefficient k, existence and uniqueness of a weak 
solution z can be proved [1]. 

In practice, for control problems u and/or g will not be prescribed, but computed as part of a 
constrained optimisation problem. 

2.2.2 Cost Junctionals 

We consider two cost functionals which, while outwardly similar, account differently for the sto- 
chastic nature of the problem. The first functional has the form: 

^ 2 /^ II l|2 

Ji{z,u,g) := -||z-z|1^2(c)oL2(o) + -||std(z)||^,^^^ 

T 2 ^2 

+ 2 ll"l'i^(-D)«>L2(n) + 2^\9\\L-^(dDN)®L^(n) ' C^) 

where 0:£)xO— >Misa prescribed target function and a, /3, 7 and 5 are positive constants. 
Typically for control problems z will be deterministic, but we permit here a more general case. 
The first term in ([7]) is a measure of the distance, between the state variable z and the prescribed 
function i, in terms of the expectation of {z — z) . The second term measures the standard 
deviation of z, which is denoted by std(z) and defined by: 



std(z) := { (z- zdV] dV] . (8) 




By increasing the value of /3 with respect to a, a greater relative contribution of the variance of z 
to J7i is implied. The final two terms in ([7]) are regularisation terms for the distributive control 
via u and the boundary control via g. 



The second considered functional has the form: 

Q'_ 2 w ||2 T2 ^2 

J2 {z,u,g) := -11^ - i||i2(c) + -||std(z)||^,^^j + ^ ll"llL2(D)»L^(n) + 2 ll5llL^(ai5„)»L2(n) > (9) 

where z is the expectation of z and z : 13 — > M is a prescribed target function. This functional is 
subtly, but significantly, different from that in ([7]). The first term in J2 measures the i^-distance 
between the expectation of z and the target z. This includes no measure of the variance of the 
actual response and is unaffected by the presence of uncertainty in z. 

Setting a = 1, /? = and 5 = Q, the functional J7i in ^ coincides with a functional presented 
in Borzi J5] , but in their work a method is constructed around a modified version of J7i . For the 
above parameters and for deterministic w, J7i coincides with the cost functional considered by Hou 
et al. [12j . 

2.2.3 Structure of the control functions 

A feature of our work is that for control problems, as opposed to inverse problems, we will consider 
control variables that are decomposed additively into unknown deterministic (to be computed) and 
known stochastic components. We will consider u to have the form 

u (x, uj) = u (x) + u' {x, U)) , (10) 

where u : Z3 — > M is deterministic and is the mean of u and u' : D x — > M is a zero-mean 
stochastic part. The goal will be to compute m, which constitutes the 'signal' sent to a control 
device. The actual controller response is w, with u' modelling the uncertainty in the controller 
response for a given instruction. The boundary control function will be decomposed analogously, 

g{x,uj)=g{x)+g'{x,uj), (11) 

where g : BD^ — !• M and g' : dDjs; x Jl — > M is a zero-mean stochastic part. 

2.3 Representation of stochastic fields 

2.3.1 Finite- dimensional noise assumption 

In representing random fields, we employ the finite-dimensional noise assumption [Tl Section 2.4], 
which states that the random fields k, zjj, g and u can be approximated using a prescribed finite 
number of random variables ^ = {^i}j^i, where L G N and ^^ : fi — >■ F^ C R. We assume that 
each random variable is independent and is characterised by a probability density function pi : 
Ti -> [0, 1]. Defining the support F — Hi^i F^ C M^, for a given y — (2/1, ■ . • , J/l) G F the joint 
probability density function of ^ is given by p = ni=i PiiVi)- The preceding assumptions enable a 
parametrisation of the problem in y in place of the random events w [S] . 

As an example, consider a finite-term expansion of the stochastic coefficient k based on L 

random variables: 

s 

i^{x,y) = ^Ki{x)(i{y) xeD,yer, (12) 

4=1 

where k^ : D -> M and C^ : F — > M. If k is represented by a truncated Karhunen-Loeve ex- 
pansion [141, then S — L + 1 with (^i = yi^i and y^ — \. If a generalised polynomial chaos 
expansion |26j is used, Q is an L-variate orthogonal polynomial of order p and S = {L + p)\ / {L\p\) . 

2.3.2 Definitions 

Given the joint probability density function p{y) of ^, with y £V and where F is the support of p, 
a space ip(F) equipped with the inner product 

{v,w)l2^t)-= vwpdy u,weLp(F) (13) 



is considered. The norm I 



\L^{D)<SLl{r] 



is then defined as: 



ll^ll 



L^D)®Ll{T) 



\z\ pdxdy 



T J D 



(14) 



Expressing two functions ui and U2 as ui{x,y) = vi{x)wi{y) and U2{x,y) = V2{x)w2{y), where 
ui, 112 € H^{D) and wi,W2 € HUT), an inner product on H^{D) (E) -ffp(r) is defined by 



{ui,U2)H^D)(>sm{r) — (v1'«2)h1(Z5) (wi,W2)ffl(r) 



(15) 



where (•, Oifiru) i^ the standard H^ inner product and (•, ■)fji is the H^ inner product weighted 
by p, analogous to (|13[). The norni||-||^i/^N„^icpN is that induced by (15). The H^ definition will 



be used in Section |4] when considering the impact of introducing extra regularisation terms into 
the cost functionals. 



2.3.3 Parametric optimal control problem 

Following from the Doob-Dynkin Lemma [Tj, the random field z can be expressed as a function 
of the given L random variables. This enables one to reformulate the stochastic optimal control 
problem corresponding to the cost functional in (l7| or (|9]) and the constraints in Q-lp]) as a 
parametric PDE-constrained optimisation problem. The parametric problem involves solving 



subject to: 



min Ji (z, u, g) or min J^2 {z,u, g) 

z^u.g z^u,g 



■ i^ix, y)W^z{x, y)) = u{x, y) 

z{x,y) = ZD{x,y) 
k{x, y)Wxz{x, y)-n = g{x, y) 



where 

Ji (z, w, g) 

and 

J2{z,u,g) 



\mD)®Ll(T) 



3td(z) 



|2 
\mD) 

, 7|| 



\mD)' 



|std(z) 



|2 



\L\D)<»Ll{r) 



\L^{D)<S,Ll{r)- 



This parametric form of the problem will be considered in the remainder. 



(16) 
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(19) 



S 2 

n\\9\\L^dDp,)(g,Ll{r) (20) 



S 2 

oll5llL2(aD„)(5L2(r) • (21) 



2.4 One-shot solution approach 



The PDE-constrained optimisation problem in ( 16 1-( 19 ) can be recast as an unconstrained opti 



misation problem by following the standard procedure of introducing Lagrange multipliers [101 123] . 
By defining a Lagrangian functional, optimality conditions can be formulated for the state, con- 
trol and adjoint variables [53] and these can be solved simultaneously. This so-called 'one-shot' 
approach yields a solution for an optimal control problem without having to apply iterative opti- 
misation routines. 



2.4-1 Lagrangian Junctionals 

Introducing adjoint variables, or Lagrange multipliers, A,x G Hq{D) (g) L^JT), we define a La- 
grangian associated with the cost functional Ji in (201 and the constraints (17)-(19) as: 



Ci{z,u,g,\x) 



r JD 



r JD 



{z- zf pdxdy+ - / z'^pdxdy-- i zpdyj dx 

A (— V2: • {kWxz) — u) pdxdy 
dz 



u pdxdy + - / g p ds dy — 
2 Jr JaDjv 



r J£) 



X 1^ 



r JdD^ 



dn 



g]pdsdy. (22) 



A Lagrange multiplier has not been introduced to impose the Dirichlet boundary condition in ( 18 1 



since this condition will be imposed by construction via the definition of the function space to 



which z belongs. For the problem associated with the cost functional J'2 in (21 1, we define the 
Lagrangian 



C2{z,u,g,X,x) ■= 



zpdy — z\ dx 



P 



+ 



u pdxdy 



r JD 



g pdsdy 



r JdD,. 



z pdx dy 



(3 



r JD 



zpdy 1 da; 



r Jd 



A (-V2: • (kVjjz) - u) pdxdy 

^{'^1 .9) pdsdy. (23) 



The control problems that we wish to solve involve finding stationary points of these La- 

grangians. The existence of Lagrange multipliers for stochastic optimal control problems is proved 

by Hou et al. [12 for Ji with a = I, /3 = and 6 = 0. The result extends to the more general 

II II 2 2 

form of Ji because of the Frechet differentiability of ||std(z)|| 2/ri-, and Ii<?|ii2/g£i )^i2/pi- 



2.4.2 Optimality system 



To find stationary points of the Lagrangians in (22 1 and (231, we consider variations with respect 



to the adjoint, state and control variables. This will lead to the first-order optimality conditions, 
which are known as the state, adjoint and optimality system of equations '10]. In the remainder 
of this section, these equations are derived by setting the directional derivative of the Lagrangians 
with respect to the adjoint, state and control variable, respectively, equal to zero. 

Following standard variational arguments, taking the directional derivative of £1 or of £2 with 
respect to x and setting this equal to zero for all variations leads to the recovery of the Neumann 



boundary condition in (19 1. Likewise, taking directional derivatives with respect to A and setting 



this equal to zero for all variations leads to the recovery of the constraint equation in (17). The 



state system of equations corresponds thus to the constraint equations (17|-(19) 



The derivative of £1 with respect to the state variable z in the direction of z* e Hq{D)(E)L'J,{T) 
reads: 

^i,z[z*]=a {z - z)z*pdxdy + (3 / zz*pdxdy-P / i zpdyj z* pdxdy 

Jr JD Jr JD Jr jd yJr J 

f f dz* 

XS/x- {K'S/xZ*)pdxdy- I / x^-i— pdsdy. (24) 

ITJD JTJdDN "■''^ 

Setting Ci^z [z*] = for all z* and following standard arguments, the following adjoint system of 



equations can be deduced: 



-Vo; • {K{x,y)\7xX{x,y)) = {a + I3)z{x,y) - az{x,y) - jS j zpdy 

K{x,y)\'xX{x,y) ■n = 
X{x,y) = 
Hx,y) = x{x,y) 

The last equation A = x is trivial. For distributive control via u the directional derivative of Ci 
with respect to u in the direction of u* E L^{D) ip(r) reads: 
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(28) 



Ci,u K] 



(7ii + \)u*p dx dy. 



r Jd 



Setting the above equal to zero for all u* implies that 

-fu{x,y) + X{x,y) ^0 xeD,yer. 



(29) 



(30) 



More specifically, considering the structure of u in (10 1, in which only the mean u is unknown, the 
optimality equation reads: 

(31) 



(32) 



ju{x) + / (7u' + X) pdy — 0. 
Since the mean of m' is zero, this optimality condition reduces to 



^u(x) 



Xpdy = X e D. 



The case of a boundary control is handled in the same fashion, but is omitted for brevity. 



The optimality conditions in ( 30 1 and ( 32 ) permit the expression of the control (u or u) as a 
function of the adjoint function A. The control can be eliminated from the first-order optimality 
equations, leaving a reduced optimality system in terms of the state and adjoint variables. In 



summary, the optimal control problem involves solving the parametric constraint problem (17) 



(19), with u or u eliminated using (30) or (32), and the parametric adjoint problem (25)-(27) 



The reduced optimality system corresponding to £2 is constructed in the same fashion as 



for £1. It consists of the parametric constraint problem in (17)-(19), with u or u eliminated using 



(30) or (32), and a parametric adjoint problem that reads: 



'Vx ■ (K(x,y)\/,^X{x,y)) = I3z{x,y) - az{x,y) + (a - /3) / zpdy 

K{x,y)'VxX{x,y) ■n = 
X{x,y) = 



x€D,yeV, (33) 

X e dDN, y e r, (34) 
X e dDo, y G r. (35) 



3 Stochastic finite element solution 



To compute approximate solutions to the optimality system derived in Section |2.4.2[ both a 
stochastic Galerkin and a stochastic collocation finite element discretisation are formulated. For 
brevity, we present the formulation for distributive control only. The boundary control case, which 
is considered in the examples section, is formulated analogously. 

For both the stochastic Galerkin and stochastic collocation methods, for the spatial domain 
we define a space Vh C Hq{D) of standard Lagrange finite element functions on a triangulation T 
of the domain D: 

Vh ■■= [vh e Hl{D) : Vh e Pk{K) V if G t} , (36) 

where K £T \s a. cell and Pk is the space of Lagrange polynomials of degree k. The space Vh is 
spanned by the basis functions {0i}j^i. 



3.1 Stochastic Galerkin finite element method 

To develop a stochastic Galerkin finite element formulation, we consider a finite dimensional space 
Yp C ip (r) for the stochastic dimension. For Yp, we adopt a polynomial basis, also known as 
a generalised polynomial chaos ^H [HI • We first consider multivariate polynomials tpg '■ T ^ M. 
generated via [16] 

L 

i}q = Yl'PqAyt), (37) 

where q = (qi, . . . jQl) € N^ is a multi-index that satisfies I?] < p and (pq. : Tg. — >■ M is a one- 
dimensional orthogonal polynomial of degree Qi. Note that p is the prescribed total polynomial 
degree and recall that L is the number of random variables in the problem. The polynomials ijjq 



are orthonormal with respect to the Lp(r)-inner product defined in (13 1. The space Yp is defined 
in terms of ipq: 

Yp := span {4,q:qeM}(lLl{r), (38) 

where A4 is the set of all multi-indices of length L that satisfy \q\ < p ioi q E Ai. It holds that 
dim (Yp) = dim (A^) = Q = { ^^) — (L +p)!/(L!p!). Consequently, there exists a bijection 

^i■.{l,...,Q}^M (39) 

that assigns a unique integer j to each multi- index /i(j) G Ai. 
A function z/jp £ V^ <E)Yp is represented as 

N 

zhp{x,y) = X! X! ^i,qMx)tpqiy)i (40) 

i=l qeM 

where Zi^q is a degree of freedom (recall that V/j is spanned by {(f>i}j^^i)- Using the above, for the 
case of the cost functional J7i with unknown stochastic u, a Galerkin formulation of the optimality 
conditions reads: find Zhp E Vh <E)Yp and Xhp G Vh<E)Yp such that 

/ / K\7xZhp-V^Whppdxdy+- / / XhpW^ppdxdy 
Jr Jd 7 Jr Jd 



gwhppdsdy V w^p €Vh<^Yp (41) 
r JdDN 



and 



KVxXhp-'^xrhppdxdy - {a + P) / Zhprhppdxdy 
ItJd JtJd 

+ /3 / i Zhppdyj i rhppdyj dx = -a / zruppdxdy "irhpEVh^Yp. (42) 



For the case that u is decomposed additively according to (10 1 and only u is to be computed 



equation (41) is replaced by: 



/ / K^^Zhp-"^xWhppdxdy+- j f / \hppdy\ f / Whppdy\dx 

= 11 u'whppdxdy+ / / gwhppdsdy. (43) 
JtJd JrJdDN 

It can be helpful to examine the structure of the matrix systems that result from the above 
finite element problems. The space Vh is spanned by N nodal basis functions, with N the number 
of spatial degrees of freedom. This enables the construction of a mass matrix M G M^^^ and 
a set of stiffness matrices Ki G M.^^^ , i = 1, . . . , S", with each corresponding to a deterministic 



diffusion coefficient Ki in (12 1. The space Yp is spanned by Q polynomials and the stochastic 
discretisation of k in ( 12 1 defines a set of matrices Ci € M'^xQ^ i — 1, . . . , S, whose elements equal 



a(j,fc) 



Ci'>P,.U)'>P,j^ik)pdy for j,k = l, 



,Q- 



(44) 



The resulting matrix formulation of the finite element problem in (41 ) and (42 1 is then given by: 



E 



a oq 

Oq c. 



>K, 






,7' 7 



(g}M 



zq 

Ai 



91 + £«! 

9q + ew'g 
—ail 

-azQ 



(45) 



where e = in the case of a stochastic control u (see (41)) and e — \ when only u is unknown 



(see (43 1). The matrix Oq e 



is a zero matrix. The diagonal matrix T S M'^xQ jg defined 



T(a, 6) :=diag[ [a a + b ... a + 6]^] 



(46) 



The vectors Zq, Xq £ M , with q = I,. . . ,Q, collect the spatial degrees of freedom of Zhp and 



Xhp, respectively, in (|4(3|; the vectors Zq,gq, u' G M.^ represent the finite element discretisation 



of /p Jjy z^qp dx dy, jpJaDjv ^i^iP ^^ ^y ''^^'^ ir Id "^'"^qP '^^ dy, respectively. 



3.2 Stochastic collocation finite element method 



In applying the stochastic collocation method [2 [3S] , we solve the optimality system at a collection 
of collocation points {y'}-_,, where y* € F. Typically, the collocation points are determined by 
constructing a sparse grid, see [31 H] for details on the point selection. Of relevance at this stage 
is that integrals of the form L {■) pdy can be approximated via 






(47) 



where Wi is an appropriate cubature weight and Q is the number of cubature points. 

From the Q realisations {z^ — 2:(a;, y') } _ , a semi-discrete global approximation of the response 
can be constructed. 



'Ax,y) ^^z{x,f)il}^{y), 



(48) 



j=i 



where the multivariate polynomials "01 are commonly interpolatory Lagrange polynomials, as de- 
fined in Babuska et al. [3]. The polynomial representation permits an exact evaluation of the 
expectation of Zp by the cubature rule (47) when the order of the cubature rule is sufficiently 
high [2F. Other moments of the response, as well its probability density function, can be easily 
extracted using a cubature rule and the expansion in ( 48 ) . 



The main computational cost of the stochastic collocation method is associated with solving 
the optimality system at each collocation point. That is, for each y*, i ~ 1,...,(3, the state 



equation (17)-(19), 

-V^ • (^K{x,f)V^z{x,f)^ = u{x,f) 

z{x,f) = ZD{x,y') 
K.{x, f)y^z{x, f)-n^ g{x, f) 



on D, 


(49) 


on ODd, 


(50) 


on dDN, 


(51) 



and the adjoint problem (25)-(28), 





E 






on a^AT, (53) 
on ODd, (54) 



are solved. For the case that the unknown control u is wholly stochastic, equation (30) is used to 
eliminate u in favour of A: 

u{x,f)^~-X{x,f). (55) 

7 



If the control function has the additive structure of (10), then 

1 ^ 
u{x) = ~-^X{x,y- 



7 



J' 



(56) 



j=i 



is used. This equation corresponds to the case where only the mean part u in ( 10 ) is to be 
computed and the stochastic part u' is prescribed. 

Remark 3.1. When moments of unknown functions, e.g., the expectation or standard deviation, 
appear in the control problem, the Q sets of equations in the stochastic collocation method become 
coupled. For the case that u = u + u' , the state equation (|49| is of the form 

(57) 



^x ■ i^K{x, f)\7ocZ{x, f)j = U{x) + U'{X, f), 



where u'{x,y^) is known. Following the usual process of applying equation (56) to eliminate the 
control function u in the constraint equation (57) in favour of the adjoint variable A, the Q deter- 



ministic problems become coupled via this term. A similar coupling is observed for /? 7^ in (52 I 



and a ^ in ( 33 ) . The advantage of decoupled systems of equations that is usually associated 



with the stochastic collocation method is therefore lost. 

Once the collocation systems are solved, the cost functional Ji can be evaluated as follows: 

Q 



Ji{z,u,g) = Vwi / Zi 



a + p 



P 



E- 



dx 
Q 



z dx 



D 



uf dx 



D 



E 



Wi I g^ ds, (58) 



with Ui = u{x, y'') and gi = g{x, y^). The scalars Wi, i — 1, . . . ,Q, represent the cubature weights 



\Q 



corresponding to the cubature points {y*}j_ij see (47). The cost functional J2 can be likewise 
evaluated. 

The coupling of the collocation systems can be visualised through a matrix formulation. Fol- 
lowing from the finite element space Vh, we define a mass matrix M g M^^^ and a set of stiffness 
matrices Ki G M^^^, i = 1, . . . ,Q, each corresponding to a K(x,y*). For a distributive control 
function u {S — in (20)), the matrix formulation of the stochastic collocation finite element 



systems in (|49[)-(|56| reads: 



7 7 



7 



Kr 



-(a + /3)/Q+/3(l^ 



-M 



Ki 



• M 



7 7 



Kr 



Z\ 



Ai 



91 + e^i 


9q + ewg 


-azQ 



(59) 
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where e = for the case that u is whoUy stochastic (see (55 1 ) and e = 1 when only u is unknown 
(see (56)). The vectors 2;^, Aj, 5^, g^, u[ € R^, with i = 1, . . . ,Q, are the finite element represen- 
tations of z{x,y^), X{x,y^), z{x,y^), g{x,y^) and u'{x,y^), respectively. The matrix Iq e M'^xQ 
is an identity matrix and the vectors 1, w G MP are defined by [1 ... 1]-^ and [wi . . . wq]"^ , 
respectively. 

In contrast with the stochastic inverse type problems in Borzi and von Winckel [7] and Zabaras 
and Ganapathysubramanian [55] to which the non-intrusivity property of the stochastic collocation 
method is preserved, deterministic simulation software cannot readily be reused for /? ^ or e 7^ 0. 

Remark 3.2. The stochastic collocation method generally requires more stochastic degrees of free- 
dom, i.e., a larger Q, than the stochastic Galerkin method in order to solve a stochastic PDE to the 
same accuracy [4 !■ Therefore, for the same accuracy, due to the coupling of the stochastic colloca- 



tion systems (see (59lj the stochastic collocation method will likely involve a greater computational 
cost compared to the stochastic Galerkin method. When presenting numerical examples, we will 
therefore only apply the stochastic collocation method to problems for which the non-intrusivity of 
the stochastic collocation method is maintained, i.e., only for the cost functional J\ with unknown 
stochastic u and P = 0. 



4 Additional regularisation of stochastic optimal control problems 

The cost functional J7i and J72 both contain a regularisation term of the form (7/2)|jw|j^2(-£))^^2(-Q-). 
This L'^-regularisation is important for the solvability of the problem, and as 7 is increased excessive 
control values are penalised [llj . However, in some applications i^-penalisation may not be 
sufhcient. We detail in this section how an additional iJ^-like regularisation can be included into 
the cost functionals and show what computational issues then follow. 

Assuming sufficient regularity of the relevant functions, we extend the functional J7i in ([7]) by 
including a iJ^-penalty on the control variables: 



Ji,m {z,u,g) 



-I 



/3| 



lL2(D),g,L2(r) 



std(z) 



|2 

\l^d) 



^11 l|2 
2\W\m{D)cSHl(T) 



5 2 

oll.9lli/i(9-D„)«>Hi(r) • (60) 



A stochastic optimal control problem involving the cost functional Ji^m and the PDE constraints 



in (17)-(19l can be solved by following the one-shot strategy described in Section 2.4 The 



state equations are given by (17|-(19) and the adjoint equations are given by (25)-(28) (see 



Section 2.4.2 1. For the distributive control u, the weak optimality condition for the Ji^h^ case 
reads: 

(61) 



T J D 



A'u*pda;d2/-h7(u,u*)^i(^)^^i(P) =0 



for all u* e H^(D) (g) Hp{T). For the additive structure of w in (10), in which only the mean u is 
unknown, the weak optimality condition reduces to 



^* + Va:U • S/xU* dx 



Xu* pdxdy =^ 



(62) 



r JD 



for all u* e H^{D), in which the zero-mean property of u' has been used. In contrast with 



the optimality conditions in (30) and (32), the optimality conditions (61) and (62) are partial 



differential equations and do not permit the straightforward expression of the control function as 
a function of the adjoint field A. As a result, no reduced optimality system is constructed and the 
optimality equations for the state, control and adjoint variables will be solved simultaneously. 
The Galerkin formulation of the optimality system is composed of three equations, which are 



the Galerkin formulation of the constraint equations (17) -(19), the Galerkin formulation of the 



adjoint equations (25)-(28), see (42), and the optimality condition in (61) or (62). Some structure 
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of the stochastic Galerkin method can be revealed from its matrix formulation. When the unknown 
control u is wholly stochastic, an algebraic system of size iNQ x iNQ results, with N the number 
of spatial degrees of freedom and Q the number of stochastic unknowns. The algebraic system 
can be shown to possess a saddle point structure of the form 



A B^ 
B 



where the matrix blocks A e m2QWx2QW ^^^^ ^ ^ ^nq^2QN ^^^^ ^^^^^ ^y. 



A = 



T(a, (i) ® M Ojvg 

OjvQ i{Iq + E)®{M + K) 



B = 



-Y.LiC^®K, Iq®M 



(63) 



(64) 



pTVxAf 



The matrix K £ 

trix T e K'5xQ jg defined in (|46|. The elements of the matrix E eMP^'^ are defined as 



is the stiffness matrix for a deterministic Laplacian operator. The ma- 



E(,3,k) := / Vj^V'mO) ■'^y'^^J.{k)P'^y j,k 



i,...g, 



(65) 



where ^^(j)(j/) denotes a multivariate, orthogonal polynomial, as defined in (37). Analytical 
expressions for computing the matrix elements in ( 65 1 are presented in Appendix | A[ 

The formulation of a stochastic collocation method follows the process detailed in Section |3.1[ 
Applying the collocation method to the optimality conditions in this section leads to Q deter- 
ministic problems which are coupled via the Q stochastic degrees of freedom. In the case of an 
unknown stochastic control function, the coupling between collocation points is due to the Vy- 



tcrnis in (61 ). When only the mean control is unknown, the discussion in Remark 3.1 is applicable 



To summarise, the stochastic collocation discretisation of the optimality system is coupled in the 
stochastic degrees of freedom when: 

• i7^(r)-control is applied; 



• the parameter /3 in ( 60 1 is non-zero; or 



• only the deterministic part of the control, u or 17, is unknown. 

In these cases, applying a stochastic collocation method is not attractive since its usual advantages 
are lost. 



5 Iterative solvers for one-shot systems 

The performance and applicability of stochastic finite element methods relies on fast and robust 
linear solvers, and this is perhaps even more so the case for stochastic one-shot methods for optimal 
control. In the context of deterministic optimal control problems, various iterative solvers have 
been designed, see for example Schoberl and Zulehner [22j, Borzi and Schulz [B] and Rees et al. 
[T5] . These solvers can readily be applied to the deterministic systems resulting from a stochastic 
collocation discretisation when there is no coupling between the collocation points. 

When a stochastic Galerkin method is applied, or in the case of coupled stochastic collocation 
systems, new iterative solution methods are required in order to efficiently solve stochastic optimal 
control problems. Such solvers typically consist of a Krylov method combined with a specially 
tailored preconditioner. This section presents two approaches for constructing preconditioners 
for the stochastic Galerkin one-shot systems. Following from Remark 3.2 solvers for coupled 
stochastic collocation systems are not considered. The presented preconditioners are applied in 
Sections |6] and [7] when computing numerical examples. 
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5.1 Mean- based preconditioner 

A straightforward and easy-to-implement preconditioner for stochastic Galerkin systems is the 



mean-based preconditioner. Apphed to the high-dimensional algebraic system in (45), the precon- 
ditioner matrix Pmcan is defined as 



n 


i' 




7 


—a 






-Pmcan := /2Q ® ^1 + ^ ®Iq®M, (66) 



where I2Q € M2Qx2Q jg ^^ identity matrix. One application of (66) requires solving Q systems, 
each of size 2N x 2iV, with N the number of spatial degrees of freedom. The solution of these 
smaller systems can be approximated by a sweep of the collective smoothing multigrid algorithm 
for deterministic PDE-constrained optimisation problems [T. 

The mean-based preconditioner performs well for problems with a low variance of n and a 
low polynomial order, as analysed for stochastic elliptic PDEs by Powell and Elman [T7]. Its 
performance deteriorates however for small penalty parameters 7 and (5, /3 ^ or e 7^ in ( |45[ ). 

Most stochastic inverse problem examples in Section [7] are solved using the mean-based pre- 
conditioner since it typically leads to the fastest solution method. The number of spatial degrees 
of freedom used in the numerical examples is small so that a direct solver for the (2A^ x 27V)- 
subsystems is most appropriate. To overcome the lack of robustness of the mean-based precon- 
ditioner, a collective smoothing multigrid preconditioner can be applied to the entire coupled 



system ( 45 1 . A suitable multigrid solver is summarised below. 



5.2 Collective smoothing multigrid 

A robust iterative solver with a convergence rate independent of the spatial and stochastic dis- 
cretisation parameters is obtained by applying a multigrid method directly to the entire system 



in (45 1 . The multigrid preconditioner is so-called 'point-based'. That is, it uses only a hierarchy of 



spatial grids. The intergrid transfer operators therefore obey a Kronecker product representation, 

hq ® Pd, (67) 

where Pd is an intergrid transfer operator based on Ki . A simultaneous update of all unknowns 
per spatial grid point is key to the multigrid performance, as discussed for stochastic elliptic PDEs 
in Rosseel and Vandewalle pS]. As a consequence, this multigrid preconditioner uses collective 
smoothing operators, e.g., a block Gauss-Seidel relaxation method. 

Most stochastic control examples in Section [6] use a collective smoothing multigrid precondi- 
tioner because of its optimal convergence properties. In contrast with the mean-based precondi- 
tioner, this multigrid method does not encounter additional difficulties when solving the coupled 



Galerkin system (45) in the case that only u is unknown (e = 1). 



6 Stochastic control examples 

A variety of numerical examples are presented to demonstrate the proposed formulation for control 
problems. When referring to control problems, we imply that the control functions are determinis- 



tic or have the additive structure in ( 10 1. Also, in the context of control we consider deterministic 
targets z only. Other scenarios are considered in Section [7] in the context of inverse problems. 

In all cases a unit square spatial domain D — (0, 1)^ is considered, and the boundary of the 
domain dD is partitioned such that dDn — {0, 1} x (0, 1) and dD^ = dD\dD£i. Unless stated 
otherwise, zero Dirichlet and Neumann conditions are applied, i.e., z^ ~ on dDjj and g — 
on dDj^. In all examples, the diffusion coefficient k is represented by a Karhunen-Loeve expansion 
based on an exponential covariancc function with a correlation length of one and variance 0.25 



(see (12)). The Karhunen-Loeve expansion is truncated after seven terms; the random variables 
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-V3,V3 



A 



i.oo.o 



Figure 1: Target function z for the control examples. 

in this expansion are assumed to be independent and uniformly distributed on 
deterministic piecewise-continuous target function is considered: 

a; e [0, 1] X (0.4, 0.6) 

1 a; e (0.1,0.9) X [0,0.4] 
z{x)={2 xe (0.2,0.8) X [0.6,1] 

10a;i X e [0, 0.1) x [0, 0.4] U [0, 0.2) x [0.6, 1] 

10 — lOxi otherwise 

This function is illustrated in Figure [T] A spatial finite element mesh of piecewise linear triangular 
elements is used. The mesh is constructed as 2'^ x 2'^ squares and each square is subdivided 
into two triangular finite element cells. This yields iV = 16 441 spatial degrees of freedom. The 
stochastic Galerkin discretisation is based on seven-dimensional Legendre polynomials of order 
two. This yields Q = 36 (when u' = in (lOl) and the algebraic system in (45) has dimension 



1.2 • 10 X 1.2 • 10 . The stochastic collocation method employs a level-two Smolyak sparse grid 
based on Gauss-Legendre collocation points with Q = 141 collocation points. 

The computer code used for the examples is built on the library DOLFIN [TB]. The complete 
computer code for all presented examples is available as part of the supporting material [20] . 



6.1 Distributive control with cost functional J'l 

First, we consider a distributive control via u only and using the cost functional J7i in ([7| with 
i5 = 0. A control function u with the additive structure in (10) is used. The goal is to determine 



an optimal mean control u, which is the deterministic input for the controller response. 



6.1.1 Perfect controller case 

The optimal control for the case u' — 0, which corresponds to the controller action being wholly 
deterministic, and /? = 0, which implies no extra control over the variance of z, is computed and 
the results are shown in Fig.[2]for the case 7 = 10~^. The computed mean of the state function, the 
variance of the state function and the control u are all shown. The same quantities are presented 
in Fig. [3] for the case 7 = lO^'^. Comparing the results in Figs. [2]and[3J as expected the larger 
7-penalty term leads to a poorer approximation of the target. The values of the cost functional, 



tracking error | 



and standard deviation std(z) 



are summarised in Table 



~||2 

The tracking error as a function of 7 is presented in Fig. [4] illustrating the deterioration in t 
quality of the computed result for increasing values of 7. 



He 



6.1.2 Imperfect controller case 

We now consider the impact of an imperfect controller by introducing a known non-zero stochastic 
term u'. The stochastic function u' is modelled by a three-term Karhunen-Loeve expansion based 
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(a) mean state 



(b) variance state 



(c) deterministic control 



Figure 2: Mean and variance of the optimal state and deterministic control variable u — u {u' = 0) 
associated with the cost functional Ji and computed with the stochastic Galerkin method with 
a=l, {3 = 5 = and 7 = 10^^. The target z is illustrated transparently in (a) for reference. 




(a) mean state 



(b) variance state 



(c) deterministic control 



Figure 3: Mean and variance of the optimal state and deterministic control variable u = u {u' = 0) 
associated with the cost functional Ji and computed with the stochastic Galerkin method with 
a = l, /3 = i5 = and 7 = 10~^. The target z is illustrated transparently in (a) for reference. 




Figure 4: Tracking error \\z — z\\j^2(jj\x^]^2(-p) as a function of the penalty parameter 7 for the 



problem considered in Section 



6.1.1 



{D)»Ll{r) 
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Xj 1.00.0 

(a) mean state 



(b) variance state 



Xj 1.00.0 

(c) mean control 



Figure 5: Mean and variance of the optimal state and mean control variable u {u' ^ 0) associated 
with the cost functional J7i and computed with the stochastic Galerkin method with a = 1, 
P = 5 = Q and 7 = 10^^. The target z is illustrated transparently in (a) for reference. 




X, "■" 1.00.0 

(a) mean state 



x^ "" 1.00.0 
(b) variance state 



X, "■" 1.00.0 
(c) mean control 



Figure 6: Mean and variance of the optimal state and mean control variable u {u' ^ 0) associated 
with the cost functional J7i and computed with the stochastic Galerkin method with a = 1, /? = 1, 
(5 = and 7 — 10^^. The target z is illustrated transparently in (a) for reference. 

on a zero-mean Gaussian field with an exponential covariance function, unit variance and a unit 
correlation length. For the examples in this section we set 7 = 10"'^. 

The computed mean and second central moment of the state z and the computed mean optimal 
control u are depicted in Fig. [5] for the case /? = 0, and in Fig. |6] for the case /3 = 1. A comparison 
between Fig. [2] which corresponds to the case u' = 0, and Fig. |5] shows that the prescribed 
uncertainty on the control has only a limited effect on the outcome of the optimal control problem 
for this example. The mean control u is visibly unchanged and the variance of the state variable 
increases slightly. This observation is also reflected in the computed values for the cost functional 
(see Table [1]). 

To provide control over the variance of the state variable, the parameter j3 in the cost func- 
tional ( pO| ) can be increased. Comparing the results in Fig. [5](/3 — 0) and Fig. [6](/3 — 1), the peak 
variance is reduced, but the correspondence between the mean state and the target is compromised. 



6.2 Distributive control vifith cost functional J2 



We now mirror the imperfect controller problem considered in Section 6.1.2 but for the cost 
functional J^2- The cost functional J7i provides a measure of the average distance between the 
state and target, whereas the cost functional J72 provides a measure of the distance between the 
mean state and the target. 

For the examples we adopt a = 1, 7 = 10"'^ and (5 = in J2. For the case /3 = 0, which 
implies no extra control over the variance of the response, the mean and the variance of the 
computed state variable z and the computed deterministic part of the control signal u are shown 
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Xj 1.00.0 

(a) mean state 



(b) variance state 



Xj 1.00.0 

(c) mean control 



Figure 7: Mean and variance of the optimal state and mean control variable u {u' ^ 0) associated 
with the cost functional S2 and computed with the stochastic Galerkin method with a = 1, 
P = 5 = Q and 7 = 10^^. The target z is illustrated transparently in (a) for reference. 




X, "■" 1.00.0 

(a) mean state 



(b) variance state 



X, "■" 1.00.0 
(c) mean control 



Figure 8: Mean and variance of the optimal state and mean control variable u {u' ^ 0) associated 
with the cost functional ^2 and computed with the stochastic Galerkin method with a — 1, /3 = I, 
6 = and 7 = 10^^. The target z is illustrated transparently in (a) for reference. 

in Fig. [7j Compared to the case presented in Fig. [5] the mean of the computed state variable 
is a better approximation of the target, while the variance of the state variable is significantly 
larger. A similar result was observed for the case u' = 0. Fig. [s] shows the computed results for 
the case /3 = 1. Increasing /3 reduces the variance of the response, but comes at the expense of 
the approximation of the target by the mean state. 

The advantage of using the cost functional J2 over J7i is that it permits a greater tuning of the 
relative importance of the mean response versus the variance of the response, via the parameters a 
and /3. In the case of the cost functional J7i, the variance of the state variable is already implicitly 
minimised by the a term as it is a norm over L^{D) x L'^{^). For the Ji case, increasing (3 only 
entails an additional contribution of the variance to the cost functional. 



6.3 Boundary control with cost functional Ji 

In practical applications it is often more plausible to control boundary values rather than the 
source term. In this section, we consider the same model problem as in Section |6.1[ but now an 
optimal boundary flux control is computed (now 7 = in the cost functionals). The right-hand 
side function in the constraint equation Q is set to m = 5. The control boundary flux is computed 
on the 'upper' (x2 = 1) and 'lower' {x2 = 0) boundaries. 
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J{z,u) 


IN 2||^2(j5)gi2(p) 


INtdW 11^(1,) 


deterministic control u = u{x), u'{x 


,y) = o 








cost functional Ji, /? = 0, 7 = 10"'^ 




2.083 X 10-' 


4.022 X 10-' 


2.562 X 10-' 


cost functional Ji, /? = 0, 7 = 10"^ 




2.911 X 10"' 


5.078 X 10"' 


1.845 X 10"' 


unknown mean control u{x), ■var(u 


)=i 








cost functional Ji, /? = 0, 7 = 10"^ 




2.956 X 10-' 


5.160 X 10-' 


1.927 X 10-' 


cost functional Ji, /? = 1, 7 = 10"^ 




3.767 X 10"' 


5.636 X 10-' 


1.367 X 10-' 


cost functional J2, /? = 0, 7 = 10"^ 




1.764 X 10-' 


2.353 X 10"' 


2.957 X 10-' 


cost functional ^2, 13 — I, -y — 10~^ 




2.956 X 10-' 


3.233 X 10"' 


1.927 X 10"' 



Table 1: Summary of the cost functional, tracking error and standard deviation of the state 
variable for the considered optimal control problems with a distributive control function. 




(a) mean state 



(b) variance state 



(c) deterministic control 



Figure 9: Mean and variance of the optimal state and deterministic control variable g — g {g' — 0) 
associated with the cost functional J'l and computed with the stochastic Galerkin method with 
a = I, -f — 0, S — 10"^ and /3 = 0. The target z is illustrated transparently in (a) for reference. 



6.3.1 Perfect controller case 

Fig. [9] presents the computed mean state, variance of the state and the boundary flux control g for 
the parameters a = l, /3 = 7 = and S = 10"'^, and with g' — 0. The corresponding values of the 
cost functional and tracking errors are given in Table [2] As can be anticipated, correspondence 
between the optimal state and target is poorer in the case of a boundary control than in the case 
of a distributive control (see Fig. Isl). The variance of the state variable in Fig. [OJca n be countered 
by increasing the parameter (3. Results for the case f3 — 1 are presented in Fig. |lO| The reduction 
in the variance is accompanied by a deterioration in the approximation of the target function. 



6.3.2 Imperfect controller case 

We now consider the impact of an imperfect controller by decomposing the boundary control g 
according to (10). The zero- mean stochastic function g' is modelled by a three-term Karhunen- 
Loeve expansion based on a Gaussian field with an exponential covariance function, a variance 
of 0.25 and unit correlation length. 

Fig. 11 presents the optimal mean boundary control g along with the first moments of the 
state variable. The values of the cost functional and tracking errors are given in Table 2] The 
mean optimal control g and the mean state appear visibly equal to the results in Fig. [9 for the 
case where a perfect control device, i.e., with g' — 0, is modelled, while the variance of z is slightly 
larger. 
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1.00.0 
(a) mean state 



(b) variance state 









— lower boundary 
""" upper boundary 


/\; 


-'/^ 







(c) deterministic control 



Figure 10: 
associated 

q;= 1, 7 = 



Mean and variance of the optimal state and deterministic control variable <7 = 5 (g' = 0) 
with the cost functional J7i and computed with the stochastic Galerkin method with 
Q, 5 = 10~^ and (3 = 1. The target z is illustrated transparently in (a) for reference. 




X, "■" 1.0 0.0 

(a) mean state 



(b) variance state 



(c) mean control 



Figure 11: Mean and variance of the optimal state and mean control variable g [g' ^ 0) associated 
with the cost functional J7i and computed with the stochastic Galerkin method with a = 1, 
/3 = 7 = and 5 = 10~^. The target z is illustrated transparently in (a) for reference. 





J{z,g) 


II "11^ 


l|std(.)||l.(^) 


deterministic control g = g{x), g'(x,y) — 
cost functional J'l, /3 = 0, 5 = 10"^ 
cost functional J'l, /3 = 1, <5 = 10"^ 


2.711 X 10"^ 
3.593 X 10"^ 


5.421 X 10"^ 
5.757 X 10"^ 


2.091 X 10"^ 
1.428 X 10"^ 


unknown mean control g{x), vsir{g')= 0.25 
cost functional Ji, P ^ 0, S = 10"^ 
cost functional Ji, /S — 1, S = 10^^ 


2.753 X 10"^ 
3.673 X 10^^ 


5.499 X 10"^ 
5.835 X 10^1 


2.168 X 10"^ 
1.506 X 10^1 



Table 2: Summary of the cost functional, tracking error and standard deviation of the state 
variable for the considered optimal control problems with a boundary control function. 
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(a) mean target 



(b) variance target 



(c) deterministic source 



Figure 12: (a)-(b) Mean and variance of the target function z for solving a stochastic inverse 
problem, (c) Deterministic source u used for constructing z as the solution of the forward problem. 



7 Stochastic inverse examples 

When the additive structure of the control function is not enforced and u is permitted to be 
stochastic and unknown, the optimal control problem effectively becomes a stochastic inverse 
problem. In this case, the stochastic properties of u are unknown, but will be computed. The 
problem is: given an observation z of a system that must obey the constraints in Q-Q, find 
the stochastic source term u that would induce the response z. Computed higher moments of u 
provide information on the uncertainty of the driving term. This formulation of the problem is not 
useful for control problems, since the properties of a control device are considered to be known and 
it is unclear how a stochastic u could be used as a control signal. The mean of u could be taken 
as the control, but this will not in general be the optimal control and computing the uncertainty 
in the system response would require an additional computation. 

Stochastic inverse problems associated with the cost functional J\ or J2 and the constraint 



equations (17)-(19) are solved using the same computational domain, boundary conditions, repre- 



sentation of K and discretisation parameters as in Section [6] For inverse problems, the collocation 
systems remain decoupled for /? = in the cost functional J\ . We present some simple examples 



using a target function z that is computed by solving the stochastic forward problem in (17)-(19l 
with the deterministic source function 



50 sin(7ra;i) cos(27ra;2). 



(69) 



Figure 12 illustrates u and the first moments of the target function that will be used. We aim to 
illustrate how the presented framework can be used for a class of inverse problems. Application 
to realistic inverse problems will require deeper investigations into handling noisy and incomplete 
data. 



7.1 Determination of the source function using cost functional J\ 

We consider the case with /3 = 0, which permits the decoupled solution of problems at collocation 
points. We therefore apply both the stochastic Galerkin and collocation methods in this section. 



1.1.1 Deterministic target 

As a first case, we consider a stochastic inverse problem with the cost functional J\ and with the 
target z taken as the mean of the forward problem. This represents an inverse problem where only 
one observation of the system response is available; some stochastic data computed in the forward 
problem has been discarded. Recall that the stochastic properties of n are known. 

The computed first moments of the state and source function, using a = l,/3 = (5 = and 
7 = 10~^, are shown in Fig. 13 The computed cost functional and tracking error are presented in 
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(a) mean state 



0.2 fJiT- 



(b) variance state 




0.2 *T. 
) 

(c) mean source 



X °'^ 1.0 0.0 



(d) variance source 



Figure 13: Mean and variance of the state and source variable for an inverse problem associated 
with the cost functional J^i and computed with the stochastic collocation method with a = 1, 
/? = (5 = and 7 = 10~^. The mean of the solution to the forward problem is used as the target. 
The deterministic source u is illustrated transparently in (c) for reference. 

Table [3j The mean of the target z and the mean of the state variable z coincide visually, whereas 
there is a considerable difference between the actual source ii and the mean of the computed 
source u. A large variance of the computed source term relative to its mean is also observed. 
This is inherent to the posed problem as limited observation data is being used. As 7 approaches 
zero, the tracking error is observed to approach zero, as shown in Fig. [14] This figure contrasts 
the tracking error in Fig. HI for a control problem, in which case the considered target cannot be 
reached as the target is not a solution of the forward problem. 



7.1.2 Stochastic target 



A stochastic target in now considered, which is the complete solution of the forward problem. 
Hence, no data has been discarded from the forward problem. The computed mean of the state 
variable, the variance of the state variable, the mean of the computed source and the variance of 
the source function are shown in Fig. 15 for the case a = l, /? = 0, J = and 7 = 10~^. The results 
were computed using the collocation method. The statistics of the state variable and the target 
visually coincide. A measure of the quality of the approximation is quantified by the tracking error 
in Table [3] As the penalty parameter 7 is decreased the computed state variable better matches 
the stochastic target. This effect is illustrated in Fig. 16 for which 7 = 10^^, and in Table [s] The 
computed source term in Fig. [16] matches the exact solution u very well in both the mean and the 
variance (which is zero for u). The stochastic Galerkin method yields similar results, as indicated 
in Table ]3] where various quantities computed with the two methods are summarised. 

For small penalty parameters, limited control over the source function is imposed, which can 
lead to a more expensive iterative solution of the one-shot systems. In the stochastic Galerkin 



case, the convergence rate of the mean-based preconditioner (see Section 5.1 ) deteriorates severely 
for small values of 7. The collective multigrid method shows robust convergence behaviour. A 
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Figure 14: Tracking error |k - ^|lL2( D),g,Lg( r) 
inverse problem considered in Section 7.1.1 



as a function of the penalty parameter 7 for the 




1.00.0 



(a) mean state 



(b) variance state 




1.0 0.0 



(c) mean source 



(d) variance source 



Figure 15: Mean and variance of the state and source variable for an inverse problem associated 
with the cost functional J\ and computed with the stochastic collocation method with a = 1, 
/3 = (5 = and 7 = 10~^. The mean and variance of the target z are illustrated in Fig. 
deterministic source u is illustrated transparently in (c) for reference. 
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1.0 0.0 



0.2 rJ,T, 



(a) mean state 



(b) variance state 




(c) mean source 



(d) variance source 



Figure 16: Mean and variance of the state and source variable for an inverse problem associated 
with the cost functional J\ and computed with the stochastic collocation method with a = 1, 
/3 = 5 = and 7 = 10~^. The mean and variance of the target z are illustrated in Fig. 12 The 
deterministic source u is illustrated transparently in (c) for reference. 
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(a) mean state 



(b) variance state 



,j^-fU, 




(c) mean source 



8 — 0.2 *T. 
3- "^ 1.0 0.0 

(d) variance source 



Figure 17: Mean and variance of the state and source variable for an inverse problem associated 
with the cost functional J72 and computed with the stochastic collocation method with a = 1, 
P = S = and 7 = 10~^. The mean and variance of the target z are illustrated in Fig. 
deterministic source u is illustrated transparently in (c) for reference. 
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The 



similar observation on the computational complexity was made by Zabaras |27| . Also, observe the 
non-smooth variance of the source function on the domain boundaries for 7 ~ 10^* in Fig. 16 'd). 



7.2 Determination of the source function using cost functional J2 



We mirror the stochastic inverse problem considered in Section [7. l.lj but now for the cost func- 
tional J^2- The ^2 formulation involves the mean of an unknown field, which means that a sto- 
chastic collocation formulation will not lead to decoupled problems when a y^ 0. Fig. [T7| shows 
the computed mean and variance of the state and source functions, computed with the stochastic 
Galerkin method. As in Fig. 13 the mean of the state variable and the mean target visually 



coincide. Since the variance of the state variable does not contribute to J72, a larger variance is 



observed than in Fig. 13 |b). The computed source function clearly differs (note the magnitudes) 
from that computed when using the cost functional J7i. The corresponding values of the cost 
functional and tracking error are summarised in Table |3] 



8 Conclusions 

A one-shot solution approach for stochastic optimal control problems with PDE constraints has 
been presented. The problem was formulated as an optimisation problem constrained by a sto- 
chastic elliptic PDE, and the framework is sufficiently general to also address a class of inverse 
problems that involve uncertainty. Statistical moments have been included in the cost func- 
tional and uncertainty in the controller response has been accounted for. To compute solutions, 
a one-shot method is combined with stochastic finite element discretisations. It is shown that the 
non-intrusivity property of the stochastic collocation method is lost when moments of the state 
variable appear in the cost functional, or when the control function is a deterministic function. We 
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J{z,u) 


II "11^ 

IF ^\\h^{D)®hl(V) 


e„ 


Stochastic Galerkin finite element solution (Q 


= 36) 






cost functional J\^ z deterministic, 7 = 10"^ 
cost functional Jx^ z stochastic, 7 = 10"^ 
cost functional J7i, z stochastic, 7 = 10"* 
cost functional J2, 7 = 10"^ 
cost functional J2, 7 = 10"* 


6.786 X 10"^ 
3.035 X 10"^ 
3.123 X 10"'' 
2.107 X 10"^ 
2.127 X 10"" 


7.225 X 10"* 
1.678 X 10"* 
1.882 X 10"*° 
3.784 X 10"® 
9.833 X 10"** 


4.368 X 10"* 
1.505 X 10"^ 
2.339 X 10"^ 
3.193 X 10"* 
3.190 X 10"* 


Stochastic collocation finite element solution (Q = 141) 


cost functional Jx, z deterministic, 7 = 10"® 
cost functional J7i, z stochastic, 7 = 10"® 
cost functional Jx, z stochastic, 7 = 10"* 


6.957 X 10"^ 
3.035 X 10"^ 
3.123 X 10"" 


7.406 X 10"* 
1.678 X 10"* 
1.882 X 10"*° 


4.556 X 10"* 
1.506 X 10"^ 
2.334 X 10"^ 



Table 3: Summary of the cost functional, tracking error of the state variable and relative error 

2 /||-||2 



U\ 



h^{D)®Ll{Y). 



\U\ 



LHD) 



of the computed source for the considered inverse problems 



with unknown source function. 



argue that the control function must contain a deterministic component in order for the problem 
to constitute a control problem, versus a stochastic inverse problem, hence the non-intrusivity 
property of the stochastic collocation method will be lost for one-shot formulations. Applying a 
stochastic Galerkin method does not impose additional difficulties compared to solving stochastic 
PDEs, hence, for the method presented in this work the stochastic Galerkin method is preferred 
over the collocation method for control problems. In the case of inverse problems, where the 
function to be found is wholly stochastic, it was shown that it is possible in some cases to preserve 
the non-intrusivity property of the stochastic collocation method. 

The formulated methods are supported by extensive numerical experiments that address both 
optimal control and inverse problems. The computed results illustrate the impact of various 
options in the formulation and the difference between the considered cost functionals. In particular, 
examples show the impact, for the considered model problem, of the different ways in which 
statistical data can be included in the cost functionals. 
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A Evaluating derivatives of stochastic functions 

Analytical expressions for inner products of the form 

j \/y^,{y)-Wyi^j{y)p<ly (70) 

are derived here based on the orthogonality properties of the generalised polynomial chaos basis 



functions ipi in ( 37 1 . Given ( 37 1 and 5 the Kronecker delta, the integral ( 70 ) can be rewritten as 

L L 






dyk dyk 



pdyk- 



(71) 



Hermite polynomials. In the case of Hermite polynomials, which are normalised with respect to 
a standard normal distribution, the derivative of ipi^ is given by 



dv>i^ 
dyk 



ikVik-i- 



(72) 
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Expression ( 71 1 then simplifies to 

/ ^^vMy) ■ V,jV'j(y)pdj/ = ^ij V jfe- 
Jr 



(73) 



fc=i 



Legendre polynomials. In the case of Legendre polynomials, normalised and scaled to a uniform 



distribution on 



-73, v^ 



it can be shown that the following relation holds: 



jfc-i 



l^-^^^f^ E v/2fe-l-2t) + l^_,_.. 



(74) 



After some calculations, expression ( 71 ) corresponds to 

L 

VyV'i(y) • "^yi>j{y)p'^y = E (^ ^ {{ik+ik) mod 2) 



A:=l 



%/2i^TTV2.7fe + 1 



min(ife, jfe) + 1 



1-2 



1 - min(ife,jfc) 



n ^n,H- (75) 
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